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Background: The trend to fabricate electrical circuits on nanoscale dimensions has led to 
impressive progress in the field of molecular electronics in the last decade. A theoretical description 
of molecular contacts as the building blocks of future devices is challenging though as it has to 
combine properties of Fermi liquids in the leads with charge and phonon degrees of freedom on the 
molecule. Apart from ab initio schemes for specific set-ups, generic models reveal characteristics 
of transport processes. Particularly appealing are descriptions based on transfer rates successfully 
used in other contexts such as mesoscopic physics and intramolecular electron transfer. However, a 
detailed analysis of this scheme in comparison with numerically exact data is elusive yet. 

Results: It turns out that a formulation in terms of transfer rates provides a quantitatively 
accurate description even in domains of parameter space where in a strict sense it is expected to 
fail, e.g. for lower temperatures. Typically, intramolecular phonons are distributed according to a 
voltage driven steady state that can only roughly be captured by a thermal distribution with an 
effective elevated temperature (heating). An extension of a master equation for the charge-phonon 
complex to include effectively the impact of off-diagonal elements of the reduced density matrix 
provides very accurate data even for stronger electron-phonon coupling. 

Conclusion: Rate descriptions and master equations offer a versatile instrument to describe and 
understand charge transfer processes through molecular junctions. They are computationally orders 
of magnitudes less expensive than elaborate numerical simulations that, however, provide exact 
data as benchmarks. Adjustable parameters obtained e.g. from ab initio calculations allow for the 
treatment of various realizations. Even though not as rigorously formulated as e.g. nonequilibrium 
Greens function methods, they are conceptually simpler, more flexible for extensions, and from a 
practical point of view provide accurate results as long as strong quantum correlations do not modify 
properties of relevant sub-units substantially. 



I. INTRODUCTION 

Electrical devices on the nanoscale have received sub- 
stantial interest in the last decade jl]. Impressive 
progress has been achieved in contacting single molecules 
or molecular aggregates with normal-conducting or even 
superconducting metallic leads [1, Hj]. The objective is to 
exploit nonlinear transport properties of molecular junc- 
tions as the elementary units for a future molecular elec- 
tronics. While the first experimental set-ups have been 
operated at room temperature, meanwhile low tempera- 
tures down to the millikelvin range, the typical regime for 
devices in mesoscopic solid state physics, are accessible 
(see e.g. This allows for detailed studies of phe- 

nomena such as inelastic charge transfer due to molecular 
vibrations voltage driven conformational changes 

of the molecular backbone [l(| , Kondo physics 11 1 , and 
Andreev reflections [6J to name but a few. 

These developments have been accompanied by efforts 
to advance theoretical approaches in order to obtain an 
understanding of generic physical processes on the one 
hand and to arrive at a tool to quantitatively describe 
and predict experimental data. For this purpose, basi- 
cally two strategies have been followed. One is based on 
ab initio schemes that have been successfully employed 
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for isolated molecular structures as e.g. density func- 
tional theory (DFT). Combining DFT with nonequilib- 
rium Greens functions (NEGF) allows to capture essen- 
tial properties for junctions with specific molecular struc- 
tures and geometries II 03 El ■ This provides insight 
in electronic formations along contacted molecules and 
gives at least qualitatively correct results for currents 
and differential conductances. However, a quantitative 
description on the level of accuracy known from conven- 
tional mesoscopic devices seems out of reach yet. Further, 
these methods are not able to capture phenomena due to 
strong correlations such as e.g. Kondo resonances. 

Thus, an alternative route, mainly inspired by solid 
state methodologies, starts with simplified models that 
are assumed to cover relevant physical features. The 
intention then is to reveal fundamental processes char- 
acteristic for molecular electronics that give a qualita- 
tive description of observations from realistic samples, 
but provide also the basis for a proper design of molec- 
ular junctions to exploit these processes. Information 
about specific molecular set-ups appears merely in form 
of parameters which offers a large amount of flexibility. 
In general, to attack the respective many body prob- 
lems, perturbative schemes have been applied, the most 
powerful of which are nonequilibrium Greens functions 
[T3. fl5|. Conceptually simpler, easier to implement, and 
often better revealing the physics are treatments in terms 
of master or rate equations. Being approximations to the 
NEGF frame in certain ranges of parameters space, they 
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sometimes lack the strictness of perturbation series, but 
have been extensively employed for mesoscopic devices 
[IH and quantitatively often provide data of at least sim- 
ilar accuracy. Roughly speaking, these schemes apply as 
long as quantum correlations between relevant sub-units 
of the full compound are sufficiently weak [TEj]. Phys- 
ically, it places charge transfer through molecular con- 
tacts in the context of inelastic charge transfer through 
ultra-small metallic contacts (dynamical Coulomb block- 
ade fPi\ ) and in the context of solvent or vibronic medi- 
ated intramolecular charge transfer (Marcus theory) fl8l — 

While rate descriptions have been developed in a va- 
riety of formulations before (2l| - [28| . the performance of 
such a framework in comparison with numerically exact 
data has not been addressed yet. The reason for that is 
simple: a numerical method which provides numerically 
exact data in most ranges of parameters space (temper- 
ature, coupling strength, etc.) has been successfully im- 
plemented only very recently in form of a diagrammatic 
Monte Carlo approach [29]. Path integral Monte Carlo 
methods have been used previously for intramolecular 
charge transfer in complex aggregates J0, [l9[ in a vari- 
ety of situations including correlated [3(| and externally 
driven transfer [3l| and, of particular relevance for the 
present work, transfer in presence of prominent phonon 
modes [32| . 

The goal of the present work is to study a simple yet 
highly non-trivial set-up, namely, a molecular contact 
with a single molecular level coupled to a prominent vi- 
bronic mode (phonon) which itself may or may not be 
embedded in a bosonic heat bath. We develop rate de- 
scriptions of various complexity, place them into the con- 
text of NEGF, and compare them with exact data. The 
essence of this study is, astonishingly enough, that rate 
theory provides quantitatively accurate results for mean 
currents over very broad ranges of parameter space, even 
in domains where they are not expected to be reliable. 



II. RESULTS AND DISCUSSION 

In Sec. Ill Al we define the model and the basic ingredi- 
ents for a perturbative treatment. A formulation which 
closely follows the P(£')-theory for dynamical Coulomb 
blockade is discussed in Sec. Ill Bl Nonequilibrium ef- 
fects in the stationary phonon distribution are analyzed 
in Sec. Ill CI based on a dynamical formulation of charge 
and phonon degrees of freedom. The presence of a sec- 
ondary bath is incorporated in Sec. Ill Dl together with 
an improved treatment of the dot-lead coupling, which 
is exact for vanishing electron-phonon interaction. The 
comparison with numerically exact data and a detailed 
discussion is given in Sec. Ill El 



A. Model 

We start with the minimal model of a molecular con- 
tact consisting of a single electronic level coupled to 
fermionic reservoirs, where a prominent internal molec- 
ular phonon mode interacting with the excess charge is 
described by a harmonic degree of freedom (cf. Fig. [T]) 
fill [33l HH . Neglecting spin degrees of freedom the total 
compound is thus described by 

H = H L / R + Ht + Hn + Hr>,Ph+ Hph 

= ^k, a c\ a c k ^+ ( T k,»c{ a d + h.c.) 

a=L,R; k a=L,R; k 

+eDd t d+ A + n^i( XQ + i Qd u) 2 (i) 

2m 2 

Here, the T& iQ denote tunnel couplings between dot level 
and reservoir a and lo = Mq ^J^jfuJ^m contains the cou- 
pling Mo between excess charge and phonon mode. An 
external voltage V across the contact is applied symmet- 
rically around the Fermi level such that €k, a = £o(fc) + Ma 
with the bare electronic dispersion relation eo(^) an d 
chemical potentials fi^ = +eV/2, = — eV/2. Below, 
this model will be further extended to include the em- 
bedding of the prominent mode into a large reservoir of 
residual molecular and/or solvent degrees of freedom act- 
ing as a heat bath. Qualitatively, since the dot occupa- 
tion d) d can only take the values q = 0, 1, the sub-unit 
Hd + Hu,pk + Hph describes a two state system coupled 
to a harmonic mode (spin-boson model |20|). Depending 
on the charge state of the dot the phonon mode is sub- 
ject to potentials V q (xo) = (muj 2 /2)(x + l q) 2 ■ Now, the 
presence of the leads acts (for finite voltages) as an ex- 
ternal driving force to alternately charging (q = 1) and 
discharging (q — 0) the dot, thus switching alternately 
between Vo and V\ for the phonon mode. The classical 
energy needed to reorganize the phonon is the so-called 
reorganization energy A = Vi(0) — Vb(0) = Mq/Swq- 
Quantum mechanically, the phonon mode may also tun- 
nel through the energy barrier located around xq — —Iq/2 
separating the minima of Vb,i- 

It is convenient to work with dressed electronic states 
on the dot and thus to apply a polaron transformation 
generating the shift Iq in the oscillator coordinate asso- 
ciated with a charge transfer process, i.e., 

U = cxp ( l P^J d \ (2) 



with momentum operator pg = iy // hmujQ/2(b Q — bo) 
where &J, 60 ar e creation and annihilation operators of 
the phonon mode, respectively. We mention in pass- 
ing that complementary to the situation here, the the- 
ory of dynamical Coulomb blockade in ultra-small metal- 
lic contacts is based on a transformation which gener- 
ates a shift in momentum (charge) rather than posi- 
tion [l7j . Now, the electron-phonon interaction is com- 
pletely absorbed in the tunnel part of the Hamiltonian, 
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FIG. 1: Single charge transfer through a molecular contact 
consisting of a single electronic level coupled to a harmonic 
phonon mode and contacted to metallic leads. Forward (no 
prime) and backward (with prime) rates are the basic ingre- 
dients for the approximate treatment, see text for details. 



thus capturing the cooperative effect of charge tunneling 
onto the dot and photon excitation in the molecule, i.e. 



H = H 



L/R 



H 



Hph 
H T 



D 



H 



Ph 



Ht with 



T k,ac{ a d exp [ ^p lo 



a=L,R; k 



+ h.C. (3) 



Single charge tunneling through the device can be 
captured formally exactly under weak conditions (e.g. 
instantaneous equilibration in the leads during charge 
transfer) within the Meir-Wingreen formulation based on 
nonequilibrium Greens functions flil [l5j . For the current 
voltage characteristics one finds 



'<"> " T / 



de- 



L^R 



x [iG > (e)-iG < (e)] 



eV 



(4) 



with energy dependent lead self-energies £ Q (e) = 
27r^ fe \Tk. a \ 2 S(e — €fc) and with the Fourier transforms of 
the time dependent Greens functions G < (t) = i(dJd(t)) 
and G > (t) = —i{d(i)cft). Upon applying the polaron 
transformation one has 



G<{t) = i/d^e~Ti polo eTi po ^ lo d(t)\ 
G>(t) = -i/e* Mt)lo d(t)d i e-% Pol °) , 



(5) 



where all expectation values are calculated with the 
full Hamiltonian ([3]). Of course, for Tfe iCt — > 0, 



the Greens functions factorize such as e.g. G < (t) 
i{<$ d(t)) £> exp[J(i)] with the phonon correlation 



iPoio p iPo(t)'o 



Ph 



(6) 



into expectation values with respect to the dot (D) and 
the phonon (Ph), respectively. Any finite tunnel coupling 
induces correlations that in analytical treatments can 
only be incorporated perturbatively. There, the proper 
approximative scheme depends on the range of param- 
eter space one considers. Generally speaking, there are 
four relevant energy scales T, L / R , Mo, k^T, and hujQ of 
the problem corresponding to three independent dimen- 
sionless parameters, e.g., 

m = - — , = LO hp , a = — . (7) 

nujQ Fkuq 

In the sequel we are interested in the low temperature 
domain 9 > 1 where thermal broadening of phonon lev- 
els is small so that discrete steps appear in the IV- 
characteristics. Qualitatively, seen from the dynamics 
of the phonon mode, two regimes can be distinguished 
according to the adiabaticity parameter E/fiwo = cr: For 
a < 1 the phonon wave packet fulfills on a given surface 
Vo or V\ multiples of oscillations before a charge trans- 
fer process happens to occur. The electron carries excess 
energy due to a finite voltage which may be absorbed by 
the phonon to reorganize to the new conformation (in 
the classical case the reorganization energy A). In the 
language of intramolecular charge transfer this scenario 
corresponds to the diabatic regime with well-defined sur- 
faces V q . In the opposite regime a > 1 charge transfer is 
fast so that the phonon may obey multiples of switchings 
between the surfaces Vq,i ■ This is the adiabatic regime. 
In this latter range the impact of the adiabaticity on the 
diabatic ground state wave functions is weak for mo < 1 
when the distance of the diabatic surfaces is small com- 
pared to the widths of the ground states. For too > 1 in 
both regimes electron transfer is accompanied by phonon 
tunneling through energy barriers separating minima of 
adiabatic or diabatic surfaces. The dynamics of the total 
compound is then determined by voltage driven collec- 
tive tunneling processes. Master equation approaches to 
be investigated below, rely on the assumption that both 
sub-units, charge degree of freedom and phonon mode, 
basically preserve their bare physical properties even in 
case of finite coupling too- Hence, since the model (H} 
can be solved exactly in the limits Too = and a = and 
following the above discussion, we expect them to cap- 
ture the essential physics quantitatively in the domain 
too < 1 and for all ratios a. We note that recently the 
strong coupling limit inclu ding the current statistics has 
been addressed as well 135, 361. 



B. Rate approach I 

The simplest perturbative approach considers the co- 
operative effect of electron tunneling and phonon exci- 
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tation in terms of Fermi's golden rule for the tunneling 
part Hj>. For this purpose one derives transition rates 
for sequential transfer according to Fig. Q] A straight- 
forward calculation for energy independent self-energies 
^l/r (wide-band limit) gives the forward rate onto the 
dot from the left lead 



T L (V,e D ) 



h 



defp 



eV 



Po(e-e D ) 



(8) 



where fp{e) is the Fermi distribution. Inelastic tunneling 
associated with energy emission / absorption of phonons 
is captured by the Fourier transform of the phonon- 
phonon correlation exp[J(t)] leading to 



^o(e) 



ifii-^'-n] 



k.i 



(9) 



with p a / e = (ttIq/2) [coth(|) =F l] denoting mean val- 
ues for single phonon absorption (a) and emission (e). 
The exponentials in the prefactor contain the dimension- 
less reorganization energy m§ = A/fuvo- Apparently, 
inelastic charge transfer includes the exchange of mul- 
tiple phonon quanta according to a Poissonian distri- 
bution. Further, one has the detailed balance relation 
Poi~ e ) — e ~ l3e Po( e )- F° r vanishing phonon-electron cou- 
pling mo — > only the elastic peak survives -Po( e ) <5(e)- 
We note again the close analogy to the P(_B)-theory for 
dynamical Coulomb blockade [17| . Moreover, golden rule 
rates for intramolecular electron transfer between donor 
and acceptor sites coupled to a single phonon mode are of 
the same structure with the notable difference, of course, 
that in this case one has discrete density of states for 
both sites [2(| [Hf . The fundamental assumption under- 
lying the golden rule treatment is that equilibration of 
the phonon mode occurs much faster than charge trans- 
fer. In the last two situations this is typically guaranteed 
by the presence of a macroscopic heat bath (secondary 
bath) strongly coupled to the prominent phonon mode. 
Here, the fermionic reservoirs in the leads impose phonon 
relaxation only due to charge transfer. Thus, for finite 
voltage the steady state is always a nonequilibrium state 
that can only roughly be described by a thermal distribu- 
tion of the bare phonon system (see below). One way to 
remedy this problem is to introduce a phonon-secondary 
bath interaction as well, see below in Sec. Ill Dl The re- 
maining transition rates easily follow due to symmetry 



T R (V,e D ) 



= ^3 
= r L (-v,-e D ). 



■T L (y,-e D ) 

T L (-V,e D ) 



(10) 



Now, summing up forward and backward events, the 
dot population follows from 



dpdo 
dt 



T to t,0 Pdot + Td 



(11) 
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FIG. 2: /^-characteristics for symmetric coupling El = E_r 
and for varying electron-phonon coupling mo at inverse tem- 
perature 6 — 25 (solid) and 6 = 10 (dashed). 



with the total rate r tot , = T L + T R + T' L + T' R and the 
rate for transfer towards the dot = + T' R obtained 
according to ((5J). Note that for vanishing electron-phonon 
coupling Mo = one has hT tot fi(M = 0) = £_[, + Sr. 
The steady state distribution pdot — > Pdot = Fd/Ft o t,0 
is approached with relaxation rate Ttot.o- For a sym- 
metric situation = E# with ejj = one shows that 
Pdot = 1/2 independent of the voltage, while asymmet- 
ric cases lead to voltage dependent stationary popula- 
tions. The steady state current is given by I(V) = 

( e /2)[(r L - ry (l- Pdat ) - (r' L - r R ) Pdot ] so that 



I(V) = e- 



rrr H -r' r r 



L x R 



tot.O 



(12) 



A transparent expression is obtained for eo = 0, namely, 
e ElSr 



I(V) 



de 



^o(e). 
(13) 



Despite its deficiencies mentioned above, the golden rule 
treatment provides already a qualitative insight in the 
transport characteristics. Typical results are shown in 
Fig. The JV-curves display the expected steps at 
eV = 2nhu>o, n 6 Z. Each time the voltage ^f- exceeds 
multiples of hu>o new transport channels open associated 
with the excitation of one additional phonon. For higher 
temperatures the steps arc smeared out by thermal fluc- 
tuations. The range of validity of this description fol- 
lows from the fact that a factorizing assumption for the 
phonon-electron correlation has been used and an instan- 
taneous equilibration of the phonon mode after a charge 
transfer, that means a < 1 and mo < 1. The latter 
constraint guarantees that conformational changes of the 
phonon distribution remain small. 

There are now three ways to go beyond this golden 
rule approximation. With respect to the phonon mode, 
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one is to explicitly account for its nonequilibrium dy- 
namics, another is to introduce a direct interaction with 
a secondary heat bath in order to induce sufficiently fast 
equilibration. With respect to the dot degree of freedom 
one can exploit the fact that for vanishing charge-phonon 
coupling the model can be solved exactly. 
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A 

a 4 - 
v 



C. Master equation for nonequilibrated phonons 
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To derive an equation of motion for the combined dy- 
namics of charge and phonon degrees of freedom, one 
starts from a Liouville-von Neumann equation for the 
full polaron transformed compound ([3]). Then, applying 
a Born-Markov type of approximation with respect to the 
tunnel coupling to the fermionic reservoirs, one arrives at 
a Redfield-type of equation for the reduced density ma- 
trix of the dot-phonon system [l5| . An additional rotat- 
ing wave approximation (RWA) separates the dynamics 
of diagonal (populations) and off-diagonal (coherences) 
elements of the reduced density. Denoting with P™ the 
probability to find q charges on the dot (here, for sin- 
gle charge transfer q = 0,1) and the phonon in its n-th 
eigenstate, one has 



dPl\t) 
dt 



E 

a=L,R; k 



\fn,k\ x 

U(EL a )K-U(-E q k OP^ 



(14) 



with vq = l,v\ = —1 and energies E%'" = hu>o(k — n) + 
v q{\J"a +££>)• The matrix elements of the phonon shift 
operator f n ^ — (n\ exp{ipolo/h)\k) read 



fn,k 



X 



( irr, \\n-k\ ( max K fc ) 

^-^i n < 

v ' y=min(ri,fc+l) 

iFi ( max(n, k) + 1, \n — k\ + 1, — m ) , (15) 



where \F\ denotes a hypergeometric function. The un- 
derlying assumptions of this formulation require weak 
dot-lead coupling a < 1 and sufficiently elevated temper- 
atures a6 < 1 for a Markov approximation to be valid. 
We will see below when comparing low temperature re- 
sults with numerically exact data that this seems to be 
only a weak constraint though. 

The calculation of the steady state distribution P™ = 

lim t _ i . 00 P" (t) reduces to a standard matrix inversion. 
One can show that for a symmetric system with (d = 
0, Si = one has Pq = P™. A typical example for the 
mean phonon number (n) = n nP™ is depicted in Fig. 
[3J The curve is well approximated by a/mo with a w 0.7. 
Apparently, (n) diverges for m — > since then Pq,P™ 
approach constants independent of the phonon number. 
Upon closer inspection one finds that excitation is more 
likely than absorption, i.e f(n, n+1) > f(n,n — l), for all 
< n < ATo(too) where iVo(mo) increases with decreasing 
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FIG. 3: Mean phonon number in nonequilibrium for eV 
3fkoo and vs. the electron-phonon coupling mo. 
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FIG. 4: Phonon number distribution in nonequilibrium for 
eV — 5hujo, mo = 0.5 and ksT/hujo = 0.1 (histogram). The 
solid line depicts a fit to a Boltzmann distribution. See text 
for details. 



mo- The opposite is true for n > Ao(too) so that in a 
steady state, depending on the voltage, the tendency is 
to have higher excited phonon states occupied for smaller 
couplings too. In particular, for strong coupling transi- 
tions n — > n + k,k > are basically blocked at quite 
small n. 

A convenient strategy to include nonequilibrium ef- 
fects in the phonon distribution, sometimes used in the 
interpretation of experimental data, is the introduction 
of an effective temperature T e ff- This way, one could 
return to the simpler modeling of the previous section. 
However, the procedure to identify PJ 1 + P" m Pa = 

exp(— /3 c fffiwo«)/[l — exp(— /3 e fffio;o)] is not reliable as Fig. 
0] reveals. While it clearly shows the general tendency of 
a substantial heating of the phonon degree of freedom 
induced by the electron transfer, the profile of a thermal 
distribution strongly differs from the actual steady state 
distribution. 

Non-equilibrated phonons leave their signatures also 
in the TV-curves as compared to equilibrated ones. The 
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pr, 



net current through the contact follows from summing 
up the transfer rates from / onto the dot according 
namely, 

1 = ^ Yl ^™< fc | 2 { ( E nk) ~ ^RfP ( E nf 

n,k 

~ [VLfp - Z R fp (E 1 * )] P?) . (16) 

Fig. [5] shows that deviations are negligible for low 
voltages in the regime around the first resonant step 
(\eV/2\ < ftwo), where at sufficiently low temperatures 
only the ground state participates so that the steady state 
distribution basically coincides with the thermal one. For 
larger voltages deviations occur with the tendency that 
for smaller couplings uiq the nonequilibrated current is 
always smaller than the equilibrated one (I n0 n < -^oq), 
while the opposite scenario (J non > I eq ) is observed for 
larger m . At sufficiently large voltages, one always has 
J non < / cq . This behavior results from the combination 
of two ingredients, namely, the phonon distributions P™ 
and the Frank- Condon overlaps \f n ,k\ 2 - To see this in 
detail, let us consider a fixed voltage. Then, on the 
one hand, for smaller mo the steady state distribution 
is broad (cf. Fig. [3]) so that due to normalization less 
weight is carried by lower lying states compared to a 
thermal distribution at low temperatures; on the other 
hand, for mo < 1 the overlaps \f n ,k\ 2 favor contributions 
from low lying states in the current (|16[) which is thus 
smaller than J eq . For increasing electron-phonon cou- 
pling mo > 1 overlaps |/ n ,fc| 2 tend to include broader 
ranges of phonon states also covered by P™ compared to 
those of low temperature thermal states. A voltage de- 
pendence arises since with increasing voltage higher lying 
phonon states participate in the dynamics supporting the 
scenario for smaller couplings. Interestingly, as already 
noted in [23| the overlaps |/ n ,fc| 2 may vanish for certain 
combinations of n, m depending on mo due to interfer- 
ences of phonon eigenfunctions localized on different di- 
abatic surfaces V q , q = 0, 1. 



D. Rate approach II 

The assumption of a thermally distributed phonon de- 
gree of freedom during the transport can be physically 
justified only if this mode interacts directly and suffi- 
ciently strongly with an additional heat bath (secondary 
bath) realized e.g. by residual molecular modes. Here 
we will generalize the formulation of Sec. Ill Bl to a situa- 
tion where the secondary bath is characterized by Gaus- 
sian fluctuations. Its corresponding modes can thus effec- 
tively be represented by a quasi- continuum of harmonic 
oscillators for which the phonon correlation function © 
can be calculated easily 
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FIG. 5: TV-characteristics for equilibrated (solid) and 
nonequilibrated (dotted) phonon distributions according to 
(113|) and (|16p . respectively. 



Here the spectral density I(uj) describes now the com- 
bined distribution of the prominent mode and its sec- 
ondary bath. It is thus proportional to the imaginary 
part of the dynamical susceptibility of a damped har- 
monic oscillator [20| . For a purely ohmic distribution of 
bath modes, one has 



= 2moCJo 



70; 



(uj 2 - LJ^) 2 + 7 



(18) 



where 7 denotes the coupling between phonon mode and 
bath. The Fourier transform of exp(J) reads at finite 
temperatures 



P 7 (e) 



-Pi 



6(e) 



-Pi 



fc! 

(M)#(0,0) 



7T 



7,e 



i\ e + nrik-hri*i 



(19) 



with frequency = u> ^ + ij/2 and /? 7i0 (f2) 



{m 2 /2£n 2 )[coth(f3hn/2) - 1] where £ = y/1 - 7 2 /4wg. 
Further, pj te (Q,) = — /9 7 , Q (— fi*) (* means complex conju- 
gation) and p 7 = K[p 7 , a +p 7i e]/2. In the above expression 
contributions from the Matsubara frequencies in Equa- 
tion (17) have been neglected since they are only relevant 
in the regime jhft ^> 2tt which is not studied here. Ap- 
parently, the coupling to the bosonic bath effectively in- 
duces a broadening of the dot levels Jvy(k+l) /2 compared 
to the purely elastic case ([9]). In the low temperature 
regime, where for equilibrated phonons absorption (re- 
lated to fc) is negligible, the widths grow proportional to I. 
The presence of the secondary bath drives the prominent 
phonon mode towards thermal equilibrium with a rate 
proportional to this broadening. Hence, if the time scale 
for thermal relaxation is sufficiently smaller than the time 
scale for charge transfer, i.e. I/77 = Sl + Sft/7 *C 1, the 
assumption of an equilibrated phonon mode is justified 
and the golden rule formulation (fl~3| can be used with 
Po(e) — > P 7 (e). However, this argument no longer ap- 
plies in the overdamped situation 7/^0 3> lj where the 



7 



phonon mode exhibits a sluggish thermalization on the 
time scale 7/wg which may easily exceed 77. 

As already mentioned above, for vanishing charge- 
phonon coupling mo = 0, the model (j3|) can be solved 
exactly to all orders in the lead-dot coupling [ll|. In 
the frame of a rate description, one observes that in 
this limit the dot population (fTTj) decays proportional to 
£l + Sfl. The golden rule version of the theory neglects 
this broadening in (|13l) since it is associated with higher 
order contributions to the current (1131) . Now, recalling 
that Po(e) reduces to a (5-function for to — > 0, this finite 
lifetime of the electronic dot level is included to all orders 
by performing the time integral in the Fourier transform 
with e -> e - i{T, L + T, R )/2 = e - iT to t(M = 0)/2 [see 
In fact, this way one reproduces the exact solu- 
tion (one electronic level coupled to leads with energy 
independent couplings), i.e., its exact spectral function. 
To be specific, let us restrict ourselves in the remain- 
der to the symmetric situation = E/j = S/2, and 
eo = 0. Then, in the presence of the phonon mode 
(too 7^ 0) the corresponding function P E ( e ) follows from 
© by replacing the S- function by i/[e + fkj(k — l) + iT,/2}. 
Again following the spirit of a rate treatment, an im- 
proved version of this result accounting for higher order 
electron-phonon correlations is obtained by using instead 
of the bare dot level width Y,/h = r to t,o(^o = 0), the 
decay rate r to t,o(-Mo 7^ 0). Equivalently, one replaces 
i/[e + hu(k + ffi/2] -> i/[e + haj(k - I) + ihT totfi /2} 
to arrive at an improved Pp(e). We note that within a 
Greens-function approach and upon approximating the 
corresponding equations of motion a similar result has 
been found in [l5|, |33| with the difference though that 
there instead of r to t,o & n imaginary part of a phonon 
state dependent self-energy E' fc ' ; appears. One can show 
that the r to t,o appearing here within a rate scheme is 
related to a thermally averaged S' fe ' ; . 

Now, an additional secondary bath can be introduced 
as above by combining (fTTJ]) with Pp, leading eventually 
to 



-P~l 



irh 



-5R 



(fe,0>(o,o) 



k\ 



n <= 



nk~n*i 



: r t , 



(20) 



The width of the electronic dot level is thus voltage de- 
pendent and approaches the bare width from below for 
large voltages, that is limy^oo r t ot,o(^) = E/h. The 
range of validity of this scheme is the following: it ap- 
plies to all ratios a /mo in the domain where the electron- 
phonon coupling is weak too < 1. For too > 1 charge 
transfer is strongly suppressed and the phonon dynamics 
still occurs on diabatic surfaces for cr/niQ <C 1 so that we 
expect the approach to cover this range as well. 



E. Comparison with numerically exact results 

A numerically exact treatment of the nonequilibrium 
dynamics of the model considered here is a formidable 
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FIG. 6: TV-characteristics according to approximate mod- 
els based on equilibrated phonons (solid) and nonequilibrated 
phonons (dashed) together with exact DQMC data (dots) 
for k&T /£ = 0.2 and without coupling to a secondary bath 
(7 = 0). All quantities are scaled with respect to the dot-lead 
coupling E. 



task. The number of formulations which allow simula- 
tions in non-perturbative ranges of parameter space is 
very limited. Among them is a recently developed di- 
agrammatic Monte Carlo approach (diagMC) based on 
a numerical evaluation of the full Dyson series which 
in contrast to numerical renormalization group (NRG) 
methods [37l | covers the full temperature range. For the 
sector of single charge transfer results have been obtained 
with and without the presence of a secondary bath inter- 
acting with the dot phonon mode. We note that com- 
putationally these simulations are very demanding as for 
each parameter set and a given voltage the stationary 
current for the IV curve needs to be extracted from the 
saturated value of the time dependent current I(t) for 
longer times. Typical simulation times are on the order 
of several days and even up to weeks depending on the 
parameter range. In contrast, rate treatments require 
minimal computational efforts and can be done within 
minutes. Here, we compare numerically exact findings 
with those gained from the various types of rate/master 
equations discussed above. 

We start with the scenario where the coupling to a sec- 
ondary bath is dropped (7 = 0) to reveal the impact of 
nonequilibrium effects in the phonon mode. The formu- 
lation for an equilibrated phonon is based on (|T3)) with Pq 
replaced by in (f2"U|) . while the steady state phonon 
distribution is obtained from the stationary solutions to 
In the latter approach the intrinsic broadening of 
the dot electronic level due the lead coupling is intro- 
duced in the following way: One first determines via (|14[) 
a steady state distribution P™ . This result is used for an 
effective self-energy contribution (total decay rate) for 
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non-equilibrated phonons, i.e., 

rVncqOO = | E l/n,*| 2 P 9 "^(^ fc ), (21) 

a— L,R; n,k,q 

where E£ k = hujQ(n — k) + /i Q . We note in passing that 
limy^oo r to t,ncq(^) = E/7i = r to t(Mo = 0). Subse- 
quently, an improved result for the steady state phonon 
distribution at a given voltage is evaluated working again 
with (|14p but using the replacement 

,neq 

J 2^ V T —) [e - hu Q (k - OP + ft 2 r t 2 ot , neq /4 • 

(22) 

Of course, for S —> the standard Fermi distribution is 
regained. The corresponding steady state phonon dis- 
tribution eventually provides the current according to 
(|T5)l using in this expression the same replacement (f2"2"j). 
The procedure relies on weak electron-phonon coupling 
mo < 1 an d requires in principle also sufficiently elevated 
temperatures. 

Results are shown in Fig. [5] together with correspond- 
ing diagMC data for various coupling strengths m . In- 
terestingly, the equilibrated model describes the exact 
data very accurately from weak up to moderate electron- 
phonon coupling too ~ 1, while deviations appear for 
stronger couplings too > 1 and voltages beyond the first 
plateau eV > 2hujQ. For mo > 1 nonequilibrium ef- 
fects are stronger and the corresponding master equation 
(HU gives a better description of higher order resonant 
steps. Moreover, as already addressed above, even in 
this low temperature domain the approximate descrip- 
tion provides quantitatively reliable results. In Fig. [7] 
the frequency of the phonon mode is fixed and only the 
electron-phonon coupling is tuned over a wider range. 
For strong coupling (here mo = 2) the equilibrated 
(nonequilibrated) model predicts a smaller (larger) cur- 
rent than the exact one in contrast to the situation for 
smaller too . This phenomenon directly results from what 
has been said above in Sec. Ill CI for stronger coupling the 
Franck-Condon overlaps favor also higher lying phonon 
states that are suppressed by a thermal distribution. 
After all, the approximate models give not only qual- 
itatively a correct picture of the exact IV curves, but 
provide even quantitatively a reasonable description in 
this low temperature domain. 

In a next step the coupling to a secondary bath is 
turned on (7 ^ 0) enforcing equilibration of the phonon 
mode, see (|20[). The expectation is that in this case de- 
partures from the equilibrated model are reduced. In 
Fig. [5] data are shown for a ratio mo = 4/5 where devia- 
tions occur at larger voltages as observed in the previous 
figures. Obviously, due to the damping of the phonon 
mode the resonant steps are smeared out with increasing 
7. However, the approximate model predicts this effect 




5 10 15 20 



V[E/e] 

FIG. 7: Same as in Fig.[6]but for fixed fk>Jo/Y, = 5 and varying 
electron-phonon coupling. 

to be more pronounced as compared to the exact data, 
particularly for stronger coupling fvy/Y, > 1, while still 
7/cjo < 1. In fact, in the limit of very large coupling 
only the k = I = contribution to (|2T)|) survives so that 
at zero temperature one arrives at 

lim /(F) = /oo-arctan [ — eF ) (23) 

1^°° 7T V" r tot,0(V / )/ 

with the current at large voltages 1^ = eS/AH and 
rtot,o(y) < where equality is approached for V —> 
00. It seems that a broadened equilibrium distribution 
of the phonon induced by the secondary bath according 
to ([20]) overestimates the broadening of individual levels. 
Since the approach is exact in the limit mo — > 0, the de- 
viations appearing in Fig. |8] are due to intimate electron- 
phonon-secondary bath correlations not captured by the 
rate approach. In the overdamped regime, i.e. 7/1^0 > 1, 
the dynamics of the phonon mode slows down and may 
become almost static on the time scale of the charge 
transfer. In this adiabatic regime an extended version 
of the master equation (fT4|) is not trivial since the con- 
ventional eigenstate representation becomes meaningless. 
One should then better switch to phase-space coordinates 
and develop a formulation based on a Fokker-Planck or 
Smoluchowski equation for the phonon. This will be the 
subject of future research. 

The essence of this comparison is that, as anticipated 
from physical arguments already in Sec. Ill Al a rate de- 
scription does indeed provide even quantitatively accu- 
rate results in the regime of weak to moderate electron- 
phonon coupling mo < 1 and for all a /mo. Deviations 
that occur for larger values of mo can partially be ex- 
plained by nonequilibrium distributions in the phonon 
distribution, where, however, the master equation ap- 
proach seems to overestimate this effect. In order to 
obtain some insight into the nature of this deficiency, 
a minimal approach consists of adding to ([T4"l) with the 
extension (|22|) a mechanism that enforces relaxation to 
thermal equilibrium with a single rate constant Tq that 
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FIG. 8: TV-characteristics in presence of a secondary heat 
bath interacting with the phonon with various coupling con- 
stants 7. Shown are approximate results (solid) using (|20[) 
and diagMC data (dots); energies are scaled with E. Other 
parameters are k^T /E = 0.2, mo = 4/5, a = 0.2. 




FIG. 9: Same as in Fig. [6] but for nonequilibrated phonons 
based on an extended master equation (solid) in comparison 
to exact diagMC data (dots). 



serves as a fitting parameter. Accordingly, the respective 
time evolution equation for Pq{t) receives an additional 
term -T [P£(t) - with the Boltzmann distribu- 

tion for the bare phonon degree of freedom Pp. Cor- 
responding results for the same parameter range as in 



Figs. |H [7] are shown in Figs. [9j [10] in comparison with 
exact diagMC data. There, the same equilibration rate 
Sro/E = 0.25 is used for all parameter sets. Astonish- 
ingly, this procedure provides an excellent agreement over 
the full voltage range. It improves results particularly in 
the range of moderate to stronger electron phonon cou- 
pling, but has only minor impact for mo < 1. The indica- 
tion is thus that electron-phonon correlations neglected 
in the original form of the master equation have effec- 
tively the tendency to support faster thermalization of 
the phonon. Indeed, preliminary results with a gener- 
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FIG. 10: Same as in Fig. [7] but for nonequilibrated phonons 
based on an extended master equation (solid) in comparison 
to exact diagMC data (dots). 



alized master equation where the coupling between di- 
agonal (populations) and off-diagonal (coherences) ele- 
ments of the reduced charge-phonon density matrix are 
retained (no RWA approximation) indicate that this cou- 
pling leads to an enhanced phonon-lead interaction and 
thus to enhanced phonon equilibration. 
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